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ABSTRACT 

A kinetic equation for Compton scattering is given that differs from the Kom- 
paneets equation in several significant ways. By using an inverse differential op- 
erator this equation allows treatment of problems for which the radiation field 
varies rapidly on the scale of the width of the Compton kernel. This inverse op- 
erator method describes, among other effects, the thermal Doppler broadening of 
spectral lines and continuum edges, and automatically incorporates the process 
of Compton heating/cooling. It is well adapted for inclusion into a numerical 
iterative solution of radiative transfer problems. The equivalent kernel of the 
new method is shown to be a positive function and with reasonable accuracy 
near the intitial frequency, unlike the Kompaneets kernel, which is singular and 
not wholly positive. It is shown that iterates of the inverse operator kernel can 
be easily calculated numerically, and a simple summation formula over these it- 
erates is derived that can be efficiently used to compute Comptonized spectra. It 
is shown that the new method can be used for initial value and other problems 
with no more numerical effort than the Kompaneets equation, and that it more 
correctly describes the solution over times comparable to the mean scattering 
time. 

Subject headings: atomic processes — radiation processes: general — radiative 
transfer 



1. Introduction 

The most fundamental approach to treating Compton scattering of photons from ther- 
mal electrons is to use a Boltzmann-like kinetic equation for the photon distribution function. 
The basic physics of the Compton process is incorporated into certain scattering functions. 
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or kernel functions, that specify the probabihty of an initial photon scattering into various 
frequency and angular ranges. This Boltzmann equation gives highly accurate results, but 
can also involve a heavy computational burden, especially when the number of scatterings 
is large. 

The Boltzmann equation can be substantially simplified when certain conditions are 
met. Consider the case where the photons and electrons arc both non-relativistic, so that 
hv/m(? ^ 1 and kT/mc^ ^ 1. Here h is Planck's constant, v is the photon frequency, m is 
the electron mass, c is the speed of light, k is Boltzmann's constant, and T is the electron 
temperature. In this case the scattering kernel is relatively narrow in frequency. In fact the 
spreading is due to the thermal Doppler width is of order Ai/ ~ i/a^/^, where. 



When the width of the scattering kernel is small compared to the scale over which the 
radiation field changes substantially, it is possible to approximate the scattering terms in 
the Boltzmann equation by a second-order differential operator acting on frequency. Such 
approximate equations are generally called Fokker-Planck equations, but for Compton scat- 
tering the term Kompaneets equation is used, in honor of its originator (Kompaneets 1957). 

In its original form the Kompaneets equation applies strictly only to homogeneous, 
isotropic radiation fields in which there is time dependence but no spatial transport. How- 
ever, since the Compton scattering process itself is not very anisotropic, the Kompaneets 
scattering term is often used for non-isotropic radiation fields, introducing some additional 
degree of approximation. In this paper, for simplicity, the transport equation is also written 
without the spatial transport terms, but it should be understood that the scattering term, 
within the isotropic approximation, is directly applicable to cases involving spatial transport 
as well. 

The Kompaneets equation has been a very useful tool for describing the process of 
Compton scattering from thermal electrons in astrophysics, when the conditions mentioned 
above are met. Second-order partial differential equations are much easier to handle numeri- 
cally than full Boltzmann equations. Within its limitations the simple Kompaneets equation 
manages to incorporate a number of important physical properties and effects, namely. 



1. Conservation of photon number. 

2. Detailed balance (correct equilibrium solution). 

3. The frequency spreading due to the thermal Doppler effect. 
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4. The frequency shift due to thermal Doppler effect (inverse Compton effect). 

5. The frequency shift due to electron recoil (Compton effect). 

6. Stimulated scattering. 

The Kompanccts equation incorporates the first two properties exactly. The remaining effects 
are accurately represented only to lowest order in the parameters hu/mc^ and a = kT/mc^. 

For all its many positive features, however, the Kompaneets equation suffers from one 
major shortcoming: it is unable to treat cases where the radiation fields varies significantly 
on the scale of Compton frequency shifts, such as can occur in the neighborhood of spectral 
lines and continuum jumps. The usual approach in such cases has been to revert to the full 
Boltzmann equation, which contains the detailed scattering kernels, or to use a Monte Carlo 
method. The simplicity of the second order differential operator is thereby lost. 

An entirely different approach to the treatment of Compton scattering was presented by 
Rybicki & Hummer (1994, hereafter RH), designed to apply to conditions typical in stellar 
atmospheres. Under these conditions the radiation fields change rapidly over the scale of 
the electron thermal frequency shift Ai/ in the neighborhood of spectral lines and continuum 
edges. The method of RH handles such rapidly varying radiation fields in a numerically 
efficient way. An essential feature of the RH method is that it expresses the radiation field 
as a differential operator acting on the emissivity, not vice versa. This departure from the 
traditional Fokker-Planck type of operator leads us to call this an inverse operator method. 

As originally formulated, the RH method incorporated only numbers 1 and 3 of the 
above physical effects, namely, photon conservation and the spreading due to the thermal 
Doppler effect. This was not a serious limitation when applied to normal stellar atmospheres, 
where these are the predominant effects. However, it seemed desirable for have a method 
that combined the advantages of the Kompaneets and RH equations, which then would be 
applicable to a much wider class of problems. 

The purpose of this paper is to present such a new kinetic equation for Compton scat- 
tering, which (like the Kompaneets equation) incorporates all the physical effects 1-6 above, 
but which also (like the RH method) has the ability to treat rapidly varying radiation fields. 
While the idea for this new inverse operator method was motivated by the RH method, the 
derivation presented here is based on a simple alteration of the Kompaneets equation, which 
is considerably simpler than the approach used in RH. The new method is only accurate to 
first order in hv/mc^ and a — kT/mc^., and it assumes that the scattering process can be 
well approximated by isotropic emission. However, since the Kompaneets equation itself has 
these same limitations, the new method cannot be judged inferior because of them. 
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Prom the standpoint of applications to stellar atmospheres, there are several advantages 
of the new kinetic equation. Like the Kompanects equation, it now incorporates the processes 
of Comptonization of the radiation field and the Compton heating/ cooling of the gas, which 
may be of importance in certain cases involving X-ray irradiation, for example. Another fea- 
ture is that it satisfies detailed balance, and thus correctly describes the approach to thermal 
equilibrium. At the same time, like the RH method, it also can handle the distortion of the 
line profiles due to Doppler broadening due to scattering on the electrons. The numerical 
implementation of the method involves only a minor modification of the RH method, which 
essentially maintains the latter's favorable timing requirements, namely proportional to the 
first power of the number of frequency points. 

For problems other than stellar atmostpheres, there are also advantages of the new 
method. For example, in studying Comptonization of X-rays, Lightman & Rybicki (1979a,b, 
1980) wrote the solution to the transfer problem as a sum of products of probability of 
scattering k times and the k-th iterated kernel function for Comptonization. As we shall 
see, the equivalent kernel function for the Kompaneets equation is not wholly positive and 
is highly singular, being a linear combination derivatives of delta functions up to order two. 
The "iterated" kernels would involve even higher derivatives, and any solution involving sums 
over such functions would be entirely unworkable. On the contrary, the kernel functions of 
the inverse operator method are quite normal functions, and their iterates can be found 
stably by numerical means, making the method of Lightman & Rybicki more practical. 

In §2 the basic properties of the Kompanects equation will be reviewed and the equations 
of the new inverse operator method will be derived. It is demonstrated that the inverse 
operator mehtod, like the Kompaneets equation, satisfies all six of the properties listed 
above. A comparison with the RH method will be given. In §3 we discuss properties of a 
rather general Boltzmann equation and show how one may define the equivalent kernels of 
any approximate version of it. This is then done for the Kompaneets and inverse operator 
methods. Numerical results are presented for the inverse operator kernel, and its properties, 
including its accuracy, are discussed. In §4 iterates of the kernel function are defined and 
discussed, and numerical examples are given. In §5 a useful formula is derived that reduces 
certain summations occuring in the formalism of Lightman & Rybicki (1979a,b, 1980) to 
numerically tractable forms. In §6 it is shown how initial value problems can be treated 
using the inverse operator method. Numerical results show the advantages of the inverse 
operator method over the Kompaneets method for short times, of the order of the mean 
scattering time. In §7 a short review is given, and some possibilities for future work are 
discussed. 
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2. The Inverse Operator Method 



For isotropic homogeneous radiation, the kinetic equation for Compton scattering can 
be written, 



Here t is the time and the photon occupation number n = n{x, t) is a function of time and 
frequency, here given in terms of the scaled frequency x — hu/kT. On the right hand side 
Ue is the electron density and cit is the Thomson cross section, which we assume is valid for 
the energies considered. Within the parenthesis the term (— n) accounts for extinction and 
the term e = e{x,t) for the scattered emission. 

The equation is completed by specifying a form for the scattering term e. In the well- 
known Kompaneets equation (Kompaneets 1957; Rybicki & Lightman 1979, §7.6) this is 
given explicitly by. 



This expression requires that the photon energy is small compared to the electron rest 
energy, hv/mc^ <S 1, and that the thermal electrons are nonrelativistic, a <^ 1. In addition, 
the radiation field n(x, t) must vary slowly with frequency x over the typical width of the 
scattering kernel. Ax ~ xa^/'^. 

The equations of the new method are found as a simple alteration of the Kompaneets 
equation. The reversion of the functional relationship (3) between n and e is easily done 
to first order in a. Since e — n io lowest order, this may be used to replace n in the term 
multiplied by ct, giving 



This implicit equation for e, combined with equation (2), represents our new formulation of 
Compton scattering. 

The Kompaneets method and the new method differ significantly in how e is determined 
from n. In the ordinary Kompaneets equation this is done by applying a differential operator 
to n, as in equation (3). However, in the new method one solves a differential equation for 
e, using equation (4). It might be characterized as applying to the radiation field the inverse 
of a certain operator related to the Kompaneets operator. For this reason we call this the 
inverse operator method. 




(2) 




(3) 




(4) 



The Kompaneets equation has the advantage that by substitution of equation (3) into 
equation (2) one obtains a single partial differential equation, which makes it much more 
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amenable to analytic treatment. In the inverse operator method one must deal with two 
separate equations, (2) and (4). However, we shall show below that for many problems this 
pair of equations is no more difficult to treat numerically than the Kompaneets equation. 

The inverse operator method based on equation (4) has very different mathematical 
properties than the Kompaneets method. For example, if the radiation field n varies rapidly 
enough, application of the Kompaneets operator can lead to negative emission terms e. This 
can be seen easily from the above equations. When the radiation field varies significantly on 
the scale of the frequency shift Ax ~ xa^^"^, the second derivative term in the Kompaneets 
operator become of order a~^n, so that the second term in equation (3) is no longer of 
order a, but is of order unity. The derivative terms, which can take either sign, are able 
to dominate the first, strictly positive, term. As we shall see below, the inverse operator 
method does not lead to such negative values of e. 

The manner of derivation of equation (4) automatically guarantees that it has the same 
order of accuracy in a as the Kompaneets equation, and thus satisfies all the properties 3-6 
above to order a. It can easily be shown that it exactly satisfies the first two properties, 
namely, photon conservation and detailed balance. 

To prove photon conservation, we multiply equation (4) by x^ and integrate over all x. 
The integrated derivative term vanishes at the endpoints, which gives. 



^0 

Multiplying equation (2) by x"^ and integrating over all x, and using equation (5), then yields 
the photon conservation result, dAf/dt = 0. 

Property 2, having the exact equilibrium solution, can be proved as follows. In equilib- 
rium, dn/dt = 0, so by equation (3) we have n{x) = e{x). Substituting this into equation 
(4) and integrating once, using the fact that e{x) — >• strongly as a; — >• oo, we find. 




(5) 



The photon density is proportional to 




(6) 




(7) 



The general solution to this is the general Bose-Einstein distribution. 



(8) 
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where the constant 7 is related to the chemical potential. This is the correct photon distri- 
bution for thermal equilibrium. 

Thus we have shown that the inverse operator equation (4) has the same general prop- 
erties 1-6 as the Kompaneets equation. Since the two equations differ only at a order of 
approximation beyond the validity of the Kompaneets equation itself, the inverse operator 
equation cannot be regarded a priori as inferior to the Kompaneets equation. As we shall 
show, the inverse operator equation is superior in some important respects. 

We conclude this section by making contact with RH. In that work the mean intensity 
J(i/) and emission integral £{1/) based on the usual specific intensity were used. These are 
related to the quantities n and e by, 

J{u) = ^n{x), E{v) = ^e(x). (9) 
c c 



In terms of these variables, (4) becomes, 

d 



-av— 



dE fhv \ c^E'' 



dp \kT J 2kTv 



+ E{v)^J{v). (10) 



It is instructive to compare this to equation (25) of RH for N — 1, which in the present 
notation is 

dE' 



d 

-aiy— 



+ E{u)^J{u). (11) 



The new inverse operator method of equation (10) is seen to be a generalization of the RH 
method for = 1 that includes extra terms involving first order derivatives of the radiation 
field. These extra terms will merely involve redefinitions of the coefficients of the second- 
order difference equations, so, from the numerical point of view, the new method is easily 
incoporated into the iterative solution scheme given in RH. 

However, now, in addition to the treatment of Doppler broadening by scattering off 
thermal electrons, all the above Compton processes 1-6, including Compton heating/cooling, 
will be properly taken into account. The new method will now be applicable to a wider 
range of stellar atmosphere problems, including Comptonization of X-rays. In cases where 
the radiation should approach thermal equilibrium, this will occur correctly. 



3. Associated Kernel Functions 



Some understanding of the relationship between the Kompaneets method and the inverse 
operator method can be gained by comparison of the equivalent kernel functions associated 
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with each method. In this section we shall define these kernel functions and give some of 
their general properties. Then we shall derive explicit forms for them corresponding to the 
two methods. 

The scattering process can be conveniently characterized by means of its kernel func- 
tion, also called scattering function or redistribution function. This function determines the 
probability that an initial photon of given frequency and direction is scattered into a final 
photon of some other frequency and direction. For the simplest case of a homogeneous, 
isotropic radiation field (the case treated in this paper), the kinetic equation can be written, 

' {[l + n{x)\K{x,x')n{x')-[l + n{x')]^K{x\x)n{x)}dx'. (12) 



UeO'TC dt J ' X^ 

The kernel function is denoted by K{x, x'). The ratio x'^/x^, which ensures the conservation 

law for photons, has its origin in certain phases space factors in the derivation (see, e.g., 
Rybicki & Lightman 1978). Stimulated scattering is taken into account by the factors 1 + n, 
which give rise to quadratic, as well as linear, terms in the radiation field. 

The integral on the right side of equation (12) must vanish when the the radiation field 
is in thermal equilibrium with the electrons, that is, when it is a Bose-Einstein distribution 
of the same temperature. The principle of detailed balance makes the stronger statement 
that, in thermal equilibrium, not only must the integral vanish, but the integrand itself must 
vanish frequency by frequency. That is, 

J2 
72 



[1 + nix)]Kix, x')nix') = [1 + nix')]%Kix', x)nix), (13) 

x^ 



for all X and x' . Substituting the Bose-Einstein form (8) into this equation leads to the 
condition for detailed balance, 

x'e^X(x, x') = a;''e^'X(x', x). (14) 

The kernel function does not depend on the radiation field n, so many of its general 
properties can be derived by considering radiation fields that are restricted in certain ways. 
In particular, it is very interesting to consider cases where n <C 1, so that the quadratic 
terms in the equation can be ignored. The kinetic equation then becomes, 

dx', (15) 



1 dn 



UgaTC dt 



K(x,x')n(x') — ^K(x' , x)n(x) 



x^ 



which is linear in the radiation field. Comparing this to equation (2) we have the general 
normalization property, 

r x'^ 

1 = / —K(x',x)dx', (16) 
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and the identification, 

4.)^jKMnWW. (17) 

The kinetic equation can also be written, 

1 dn 



UeCTTC dt 



-n+ / K{x,x')n{x')dx'. (18) 



The most significant distinction between the Kompaneets method and the new method 
presented here is found in the kernel functions associated with each. The kernel functions 
associated with each method are not immediately obvious, but there is a simple trick for 
determining them based on equation (17), which applies to the case of no stimulated emis- 
sion. After first changing the dummy variable of integration to x" (say), we note that by 
substituting the special radiation field n(x) = 5{x — x') into the integral, we obtain the kernel 
function e{x) — K{x,x'), where x' is considered as a fixed parameter. 

We shall now discuss these kernel functions, first for the Kompaneets equation, then for 
the new method. 



3.1. The Kompaneets kernel 

Using the trick given at the end of the previous section, the kernel function Ki<i{x,x') 
corresponding to the Kompaneets equation can be found by substituting n{x) — S{x — x') 
into equation (3) (without the stimulated term), which gives, 

Xk(^, x') = 5(x - x') + 4|- [x'{S'{x - x') + 5(x - x')}] . (19) 

x'^ ox 

Performing the differentiation, we can also write the Kompaneets kernel in the form, 

Xk(x, x') = (1 + Aax)5{x - x') + a{x'^ + Ax)5'{x - x') + ax'^5"{x - x'). (20) 

This shows that the Kompaneets kernel is a linear combination of generalized functions, 
namely delta functions and derivatives of delta functions up to order 2, which keep it con- 
centrated near the point x = x'. In terms of a limiting process used to define the delta 
function, such derivatives of delta functions take both positive and negative signs, and so in 
a real practical sense the Kompaneets kernel is not a completely positive function. This can 
also be verified by evaluating equation (17) for the Kompaneets kernel with a positive, but 
rapidly varying, function n, say a sufficiently narrow Gaussian, which will produce a result 
e that can take negative values. 



An important property of the kernel is its set of frequency moments, which characterize 
the frequency of the scattered photon relative to the frequency x of the incident photon. 
These are defined as, 

{{x - x'Y) - Ji^- ^y^K{x, x') dx, (21) 

where s = 0, 1, 2, . . .. The factor x^ jx'"^ is included to convert the kernel K into the appro- 
priate kernel for describing photon numbers, rather than occupation number. Substitution 
of the Kompaneets form for the kernel (19) and integration by parts, one obtains the general 
result, 

((x - x'Y) = 5,0 + OL{^x' - x'^)5,i + 2ax''^?>s2- (22) 

In particular, we have. 



{{x-x'f) 


= 1, 


(23) 


{{x-xr) 


= a{4x'-x'^), 


(24) 


{(x-xr) 


= 2q;x'^ 


(25) 


{(x - xr) 


= 0, for s > 3. 


(26) 



The first of these is a simple restatement of the photon conservation law for scattering. The 
second and third give the results for the mean shift and variance of the emitted frequency 
relative to the incident frequency, respectively. These values for the moments are, of course, 
well known, and were used in the derviation of the Kompaneets equation. The above deriva- 
tion simply verifies that the Kompaneets kernel we have obtained is consistent with these 
values. 

The final result (26) states that all moments of order 3 or greater vanish. This is a 
consequence of the derivation of the Kompaneets equation, in which a Taylor series in the 
emitted frequency is only carried out up to second order, implicitly assuming that higher 
order moments would vanish. 



3.2. The Kernel for the Inverse Operator Method 

We have seen that the Kompaneeets equation has an equivalent kernel that is a gener- 
alized function, containing up to second derivatives of delta functions. We now want to find 
the equivalent kernel for the inverse operator method for treating Comptonization. First 
of all the equation for the method will be written without the stimulated scattering terms, 
since, as we have seen, the kernel function can more easily be found for this purely linear 
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equation. Dropping the term in equation (4) gives, 



n 



e — a 



]_d_ 
dx 



X 



OX 



(27) 



The kernel function determined by this equation will be denoted by K{x, x'). Applying 
the same trick as before, we set n{x) — S{x — x') and e{x) — K{x,x') in equation (27), which 
yields. 



1 d 

S{x - x') = K - a— — 
x^ ox 



X 



ox 



(28) 



That is, K — K{x, x') is the Greeen's function of the differential operator on the right hand 
side of equation (27). 

It is actually possible to solve this equation analytically in terms of Bessel functions, but 
we omit details of this derivation. We only quote one analytic result for very small values of 
a. Defining, 



9 1 

( 7 + - 

4 a 



1/2 



(29) 



the kernel is asymptotically. 



2an \x 



x^ 



(30) 



Here a;> and a;< denote the larger and smaller of x and x', respectively. Because n 



a' 



-1/2 



^ 1, the last factor dominates the behavior of the kernel, making it very sharply 
concentrated about x = x'. By virtue of the slowly varying factors (x' /x^Y^'^ exp[(a;' — x)/2], 
this approximate form of the kernel still obeys the detailed balance relation (14). 

When the kernel is sufficiently sharp, one may set x = x' in the slowly varying factors. 
Taking into account the differing notations, the kernel (29) then reduces to the form of the 
kernel given in equation (17) of RH for = 1, showing again the relation of the inverse 
operator method to RH. 

It is of great practical interest here, as it was in RH, that the scattering function e, or 
the kernel function A', is the solution of a second-order differential equation, (27) or (28), 
respectively, since these arc easily solved numerically. The numerical solution begins with 
the introduction of a frequency grid of A^; points. Using an obvious matrix notation over 
the frequency space, the Kompaneets relation (3) becomes. 



e = (1 + Q;T)n, 



(31) 
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where 1 is the unit matrix and T is a tridiagonal matrix formed from coefficients of the 
discretization of the second-order differential operator. 

Similarly, the inverse operator relation (4) becomes, 

n = (1 - aT)e, (32) 

which can be solve for e in terms of n, 

e=(l-aT)-^n. (33) 

Comparison of equations (31) and (33) should again make clear why we call this an "inverse 
operator" method. The numerical solution of a tridiagonal system of equations is quite 
rapid, with timing of order A^^:- This is the same order as the multiplication of a tridiagonal 
matrix times a vector, so that the numerical burden associated with the new operator for 
many applications is of the same order as that of the Kompaneets operator. 

The inverse operator kernel is similarly found by discretization of equation (28), giving, 

D = (1 - aT)K. (34) 

Here D is a discrete matrix version of the delta function S(x — Xq). This is a diagonal matrix 
with diagonal elements that depend on the quadrature scheme to evaluate integrations over 
X (for the crudest scheme these are equal to the inverse of the local frequency differences). 
Equation (34) can be solved column by column in order operations to give the complete 
inverse, 

K = (1 - q;T)-^D. (35) 

If the kernel is only wanted for a few values of Xq, the number of operations is correspondingly 
smaller. 

In figure 1 results for the inverse operator kernel function K{x, xq), obtained numerically 
in the above way, are plotted as the solid curves for a = 10^'^ and for values of xq = 0.1, 1, 
and 10. The normalizing factor (x^/xg) is chosen to give equal areas under all the curves. 
Results for the asymptotic results (29) are also plotted as solid curves, but these are so close 
as to be indistinguishable. As expected, the kernel is skewed towards larger (smaller) values 
of x for xq <^ 1 {xo ^ 1), expressing the tendency for the radiation to approach thermal 
equilibrium. This can be seen most clearly for the Xq — 10, less so for xq — 0.1. 

In order to get some idea of the absolute accuracy of the inverse operator results, we 
compare them to an expansion of the true kernel to third order in a given by equation 
(19) of Sazonov & Sunyaev (2000). These results are plotted as dashed curves in figure 1. 
These show that the relative errors in the inverse operator are less than 20% in a central 
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region of the kernel about xq, which extends out to values of frequency x where the kernel 
itself falls to a few percent of the central value. Beyond this central region the relative errors 
increase rapidly, since the true kernel decays faster than any power law. However, we remark 
that the absolute errors outside the central region are actually quite small. More extensive 
comparisons show that the above errors are characteristic of all the cases for which a < 10~^ 
and x/xo < 10. 

One can imagine cases where the large relative errors outside the central region could be 
a problem, for instance, for accurately describing the wings of a Compton-broadencd line far 
away from line center. However, it is not common for the kernel to appear in such pure form, 
and more typically the small absolute errors in the kernel will cause lesser problems when 
integrations over frequency are done; recall that the kernel does have the proper frequency 
moments up to second order. 

The errors in the inverse operator kernel are larger when a > 10~^, where additional 
relativistic corrections, omitted in the Kompaneets equation, begin to be important. The 
errors are also larger when xq > o;"^/^, since the dispersion associated with Compton recoil is 
not treated adequately by the Kompaneets equation, as first discussed by Ross et al. (1978). 
However, the Kompaneets equation suffers from the same problems in these regimes. 

The inverse operator kernel shares some important properties with the exact kernel: It 
is positive and has roughly the right shape in the central regions, including the characteristic 
discontinuity of slope for x — Xq. Given that no physics beyond the Kompaneets equation 
has been assumed, even the modest accuracy of the inverse operator kernels seems most 
fortuitous, especially as the Kompaneets kernel itself is singular and not wholly positive. 



A concept that often proves useful in practice is that of iterated kernels. The physical 
interpretation of equation (17) is that the kernel function K{x,x') describes how photons 
of initial frequency x' are redistributed to other frequencies x after one scattering. The 
distribution after two scatterings is described by a function K2{x,x'), clearly given by. 



and, more generally, after k > 2 scatterings, the distribution is described by a function 
Kk{x,x'), which is defined by the recurrence. 



4. 



Iterated kernels 




(36) 




(37) 
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For convenience, the definition of iterated kernels may be extended to A; = 1 and k — hy 
the formulas, 

Ki{x,x') ^ K{x,x'), Ko{x,x') ^5{x-x'). (38) 

It follows that for any nonnegative k and I, 

Kk+i{x, x') = j Kk{x, x")Ki{x", x') dx". (39) 

These formulas may, in principle, be applied to the Kompaneets kernel (20). However, 
without working out the details, it is obvious that the A;-th iterated kernel is highly singular, 
involving derivatives up to order 2k of delta functions, and is not wholly positive. Con- 
sequently, even for moderate values of k and relatively smooth functions, finding iterated 
Kompaneets kernels by the above recurrence relations will be impractical numerically. 

For the inverse operator method, the iterated kernels, like the kernel itself, will be ordi- 
nary positive functions, which can be easily computed numerically. Using the discretization 
introduced above, wc note that the interated kernel for the inverse operator method can be 
expressed as the recurrence relation 

Kfc = (1 - Q;T)-iKfe_i. (40) 

The application of the inverse of a tridiagonal matrix involves of order A^^, opperations, so 
the recurrence relation can be apphed using of order Nl operations. This is an order than 
the naive A^^ operations if full matrices were involved. If one only wants the A;-th scattered 
emission, 

ek{x) ^ J Kk{x,x')n{x')dx', (41) 

then the iterative step, 

ek{x) = J K{x,x')ek-i{x')dx', (42) 
can be formulated, in matrix language, as, 

efe = Kefc_i = (l-aT)-^efc_i, (43) 

which can be computed in order operations. 

The iterates of the kernels in figure 1 up to A; = 5 have been computed in this fashion 
and are shown in figure 2. The discontinuity in slope of the kernel itself disappears after 
one iteration, and the higher iterates become broader and quite smooth. The bias towards 
approach to thermal equilibrium can be seen in the tendency for the higher iterates to move 
toward the central Wien frequency of x ~ 3. 
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Repeated application of equations (40) and (43) lead to the explicit forms, 

Kfc = (l-aT)-^D, (44) 

and 

efc = (l-aT)-^eo, (45) 
in terms of powers of the matrix [1 — aT)~^ . 

5. A Summation Formula 

Equation (45) is the basis of a formula we have found very useful for the computation of 
Comptonized spectra using the formalism of Lightman & Rybicki (1979a,b, 1980). Putting 
that formalism into the present notation, a spectrum s{x) is expressed as the sum, 

oo 

s{x) = J^Pfeefe(a;), (46) 

k=0 

where pk is the probability of a photon in the spectrum having scattered k times before 
escaping. Again, in matrix form. 



X^Pfeefc. (47) 



k=0 

It turns out that pk can often be accurately expressed as the linear combinations of 
exponential terms of the form Az''^, where A and z > 1 are constants. We shall derive a 
formula for the spectrum associated with just one term, since the large-A; behavior of Pk is 
often dominated by one term; examples of this can be found in Nishimura et al. (1986). In 
any case, by linear superposition, we can use the formula to find the spectrum due to a sum 
of such terms. 

Thus we consider the exponential expression for p^, 

Pk = Az-'' = {l-z-')z-', (48) 

where for convenience we take ^4 = (1 — z~^) so that the pk are normalized. 



fe=0 



The summation (47) for this choice of Pk and with equation (45) is, 

oo 

s = Aj2 [^"'(1 - «T)-i] ' efc = A [1 - z-\l - aT)-'] e^, (50) 



k=0 
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summing the geometric series. A slight rearrangement yields the desired formula, 

s = ^eo + ^^"^(l-aT-;2-^)-^eo. (51) 

The second term may be found as the solution of a tridiagonal system in order A^-^ operations, 
so the computation is very fast. 

Examples of spectra computed using our summation formula are given in figure 3 for 
a source eo(a:) of soft photons, emitted in a Wien law of temperature Tq — 10^ K, which 
are Comptonized by scattering in a hot medium of temperature T — 10^ K, making a — 
1.68 X 10~^. The choice of normalizing factors here for s{x) makes the plot analogous to 
a loguFi, vs. logz/ plot, where F^, is the photon flux. The values of z are conveniently 
characterized by the parameter, 

y^ = a/lii.z, (52) 

which is a kind of overall "^/-parameter" for the scattering process. To see this, we examine 
the limit where z is close to unity, so that Inz ^ z — 1. But the mean number of scatterings 
is (A^) = ^j^. kpk = {z — 1)~^ from equation (48), so that in the limit we have, 

y,r^a{N), (53) 

which is a typical definition for a ^/-parameter; see, e.g., Rybicki & Lightman (1979, Eq. 
7.41a). In figure 3, as the values of range from 0.03 to 10, the spectrum changes from 
being nearly a Wien law at the initial 10^ K temperature to being nearly a Wien law at 
the temperature 10^ K of the Comptonizing medium. For values of frequency between these 
two Wien laws, one sees an approximate power law, which is originially steep, but which 
eventually becomes flat. 



6. Initial value problems 

In an initial value problem the radiation field is specified at an initial time, say t = 0, 
and the problem is to find its subsequent behavior as a function of frequency and time, 
that is, given n{x, 0), find n{x,t). In this section we shall describe the results of numerical 
solutions of initial value problems for the Kompaneets and the inverse operator methods for 
the spatially homogeneous case. This will help clarify some of the distinctions between the 
two methods. 

For simplicity the stimulated scattering process will be neglected, so that the equations 
are linear. The prototypical initial value problem is to start with an initial radiation field 
concentrated at a single frequency xq, so that n{x, 0) = 5{x — xq). The discretized vector 



-17- 



version of this initial function will be denoted no. Any initial value problem can then be 
solved by linear superposition of such solutions. 

It is convenient to replace t by the Compton y -parameter, 

y^ar^ UeCTTCt, (54) 

\mc^ J 

which measures time in units of the time for significant Comptonization to take place. Com- 
bining equations (2) and (3) and writing the result in matrix form, we then have, 

|^T„. (55) 

This equation can be solved by the well-known Crank- Nicholson method. Introducing a 
discretization in y (time), indicated by a subscript j, we write. 

That is, we use a differencing scheme that is semi-implicit in time. This can be rearranged 
into the form, 

(1 - AoT)n,+i = (1 + AoT)n,-, (57) 

where, 

Ao = Ay/2. (58) 

For the initial value problem we are given the radiation field no at 7/ = 0, and we con- 
struct the n's at succeeding times by using (58) as a recurrence relation. Given the solution 
n^ we evaluate the right hand side (this involves multiplying by a tridiagonal matrix) then 
solving the resulting equation for n^+i (this involves solving a linear system with tridiagonal 
matrix). Both operations involving tridiagonal matrices requires only N^. operations, so they 
are very rapid. 

Now let us turn to the initial value problem using the inverse operator method. Using 
the same discretization, equations (2) and (4) become, 

^ = a-^ \-n + (1 - aT)-^nl = (1 - aTY^Tn. (59) 
dy 

This can be differenced using the Crank-Nicholson method, as before, so that, 

^^^^ = \{1 - aT)-^Tn, + 1(1 - aT)-^Tn,+,. (60) 
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After multiplication by (1 




(61) 



where, 



A. 



Ay/2 + a, 



A_ = Ay/2 - a. 



(62) 



Equation (61) is a trivial modification of (57), the only change being the replacement of the 
two appearances of Aq with A+ and A_ on the opposite sides of the equation. Thus the 
numerical method developed for the Kompaneets equation can be trivially adapted to solve 
the initial value for the inverse operator method, and it will be equally fast. 

One distinction between the two methods is now apparent. For the Kompaneets equa- 
tion, the time variable t and a can be combined into a single variable, the y-parameter. This 
implies that all initial value problems starting from the same initial radiation field will have 
exactly the same evolution for different values of a, except for a simple rescaling of the time 
variable. For the inverse operator method, even using the y-parameter does not completely 
eliminate the dependence on a, as can be seen from equation (59), or from the fact that a 
appears in the definitions of A_|_ and A_. This implies that the inverse operator method 
describes the solution not only on the time for significant Comptonization, but also on the 
typically much shorter mean free time for scattering. 

An example of the numerical solution of initial value problems is shown in Figure 4. 
Here a = 10~^, and the initial radiation field is a delta-function at the frequency Xq = 10~^. 
The values of the y-parameter range from to 4 in steps of 0.1. For the larger values of y, 
the solution is seen to be converging to the Wien law, (l/2)a;^ exp(— x), shown as the dotted 
curve. Both the Kompaneets and the inverse operator results are plotted in figure 4, but for 
this set of y values the two solutions cannot be distinguished to plotting accuracy. 

The discussion above suggests that the distinctions between the two methods will be- 
come manifest for times comparable to the mean free scattering time, t < {ueCrTc)'^, that 
is, when r < 1 or y < a. 

Figure 5 demonstrates that the behaviors of the two methods at relatively short times 
are indeed very different. The upper set of curves show the inverse operator results and the 
lower set the Kompaneets results. To avoid confusing overlapping curves, the times have 
been divided into two sets, r = 0.25 to 1.5 in steps of 0.25 for the left panels (a), and r = 1.5 
to 15 in steps of 1.5 for the right panels (b). The initial radiation field n{x, 0) = 5{x — Xq) 
with Xq — 10~^ has not been plotted here, but would have appeared as a vertical line in all 
plots. 



For the inverse operator method, one component of the solution is clearly a remnant of 
the localized initial delta function. Immediately outside this, there is an extended component 
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that increases with time until r ~ 1.5, then falls. The extended component originally has a 
noticeable discontinuity of slope at x — Xq, but eventually becomes smoother as it broadens. 

This behavior of the inverse operator results is easily understood analytically. One 
shows by direct substitution that the solution to the time dependent equation (18) with 
n{x, 0) — S{x — xo) can be written in terms of iterated kernels, 

k=0 

For moderate times r < 1, we expect only a few terms of this series to contribute, so that 

n{x, t) = e~'^5{x — xq) + Te~'^K{x, xq) + . . . . (64) 



Thus, for the inverse operator method, we expect an early time behavior of a decaying 
delta function plus a rising contribution of the noniterated kernel, peaking at r = 1, then 
decaying away. The presence of the delta function expresses the fact that after a time r, 
some fraction (in fact, e""^) of the photons have not yet scattered; in principle, this delta 
function never goes away, but in the numerical solution it eventually becomes negligible. 
The presence of the noniterated kernel can be recognized by the discontinuity in slope at 
Xq, but at later times this feature decays and disappears as further terms involving iterated 
kernels begin to dominate. This is precisely what is seen in the numerical solutions for the 
inverse operator method in the upper panels of figure 5. 

Equation (63) is not a useful representation of the solution of the Kompaneets equation, 
due to the singular nature of the iterated kernels, so the preceding discussion is not applicable 
to it. The numerical solution of the Kompaneets shown in the lower panels of figure 5 never 
exhibits a trace of a delta function for r > 0, but rather shows quasi- Gaussian forms, which 
continually broaden and shift with time. For the largest time (r = 15) the results are 
virtually indistinguishable from those of the inverse operator method. 

The coefficient of the iterated kernel in the representation (63), that is, T^e~'^/k\, has 
the same form as the Poisson probability distribution, and for large r (or large k) is sharply 
peaked about k = t. This suggests the large-/c asymptotic approximation to the iterated 
kernels, 

Kk{x, Xq) ~ n{x, T ^k), (65) 

where n{x, r) is the initial value problem with n{x, 0) = 6{x — xq). Indeed, this is essentially 
the basis of the method of Sunyaev & Titarchuk (1980), who used it to determine Compton 
spectra. 



20 



For small r the mixing of the different values of k in equation (63) is significant, and the 
approximation (65) is not so good. Nonetheless, Sunyaev (1980) suggested using equation 
(65) for A; = 1 to define the approximate kernel. 



which he called the "Kompaneets kernel." This differs from our "Kompaneets kernel" Kk, 
which is the true kernel of the Kompaneets equation, as given in equation (20). The kernel 
Kssix,XQ) does have the desired property that the k-th iterate of it is exactly equal to 
n{x, T = /c), so it does give asymptotically accurate iterates through equation (65). However, 
we find that Kss is not as accurate as the inverse operator kernel in the central regions, and it 
does not have the characteristic discontinuity of slope at x — xq. In addition, Kss is actually 
harder to compute than the inverse operator kernel, since it is defined by an integration of 
the time-dependent Kompaneets equation over a finite interval of time. 



We have demonstrated a simple modification of the Kompaneets method that maintains 
most of its desirable properties, but which has much better behavior at very short times, 
of order of the mean free scattering time. This inverse operator method has an equivalent 
kernel function that is positive and nonsingular, unlike that of the Kompaneets equation. 
The kernel function is roughly of the right shape in the central regions with an accuracy of 
about 20%. We have shown that for many applications the inverse operator method requires 
no more numerical effort than the Kompaneets equation. 

We have shown how the iterated kernels of the inverse operator method can be efficiently 
computed numerically. A summation formula involving iterated kernels is given that reduces 
the effort in computing Comptonized spectra using the method of Lightman & Rybicki 
(1979a,b, 1980). We have shown that initial value problems can be solved using the inverse 
operator method with no more effort than the Kompaneets equation, but with noticeable 
improvement in the solution at times comparable to the mean scattering time. 

In this paper there has been no attempt to advance beyond the physics of the Kom- 
paneets equation. The main goal has been simply to point out the special advantages that 
occur when one reverses the role of the second order differential operator that connects the 
radiation field with the emission coefficient. However, in future work it would be desirable 
to overcome some of the limitations of the inverse operator method, as others have done for 
the Kompaneets equation. There does not seem to be any reason to expect any difficulties 
in extending the present theory to anisotropic scattering or polarization. 




(66) 



7. Summciry and Discussion 
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It would be desirable to find ways to increase the accuracy of the inverse operator 
method. Obvious improvements would come with using higher order expansions of the 
coefficients in the parameters a and hu/mc^, but this will not solve all of the accuracy 
issues. 

More improvement might be gained by expressing the emission coefficient as a sum of 
N > 1 terms, as in RH. Alternatively, and more in keeping with the present approach, one 
could include higher order moments and derivatives in a Kramers- Moyal expansion (see, e.g., 
Riskin 1984) of the Boltzmann equation, and then use this to form new, higher order, inverse 
operators. 

The present work has concentrated on one particular Fokker-Planck equation, the Kom- 
paneets equation. However, the idea of using inverse operators clearly could be generalized 
and applied to other Fokker-Planck equations. Since each physical situation has its own 
special properties, it is not possible to predict whether an inverse operator approach would 
provide significant advantages in any particular case, but it might be worth investigating. 
An intriguing possibility is that higher order Kramers-Moyal expansions might have better 
properties when modified by the use of inverse operators, which perhaps might avoid some 
of the difficulties known to exist for the traditional expansions beyond the Fokker-Planck 
approximation (Pawula 1967). 

The author wishes to thank A. Loeb and R. Narayan for helpful discussions and com- 
ments. The paper was also improved through the valuable suggestions of S.Y. Sazonov, D.G. 
Hummer, and especially the referee, I. Hubeny. 
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Fig. 1. — Kernel functions for xq = 0.1, 1, and 10, and for a = 10"^. The solid lines 
are the inverse operator kernels. The dashed curves are accurate approximations based on 
expansions of Sazonov & Sunyaev (2000). 
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Fig. 2. — Iterates Kk{x,Xo) of the inverse operator kernels of figure 1 for A; = 1, . . ., 5. The 
narrowest curves are the noniterated kernels (fc = 1). The iterates become broader as k 
increases. 
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log X 

Fig. 3. — Comptonized spectra from a soft initial Wien source (Tq = 10"^ K) scattering in a 
hot Comptonizing medium (T = 10^ K) with simple model for the escape probability (see 
text). The solid curves, progressing from left to right, are for values of|/* = 0.03, 0.1, 0.3, 1, 
3, and 10. The Wien curves for 10^ K and 10^ K are shown as dotted on the left and right, 
respectively. 
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Fig. 4. — Solutions (solid) for the initial value problem with an initial delta-function at 
xq — 10~^ (vertical hne) and for y = 0, 4(0.1). For large values of y the solutions ap- 
proach the Wien law (dotted) . The Kompaneets and inverse operator solutions are virtually 
indistinguishable in this plot. 
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Fig. 5. — Short time behavior of the initial value problem with a = 10^^ and an initial delta- 
function at Xq = 10~^ for the inverse operator method (upper panels) and the Kompaneets 
equation (lower panels). Values of r are (a): 0.25,1.5(0.25) for the left panels and (b): 
1.5,15(1.5) for the right panels. Note the linear vertical scale here. 



